#' ---
#' title: "Exporting Descriptive Distributions"
#' author: "Gento Kato & Fan Lu"
#' date: "June 12, 2023"
#' ---
#' 
#' 
#' # Preparation 
#' 

## Clean Up Space
rm(list=ls())

## Set Working Directory (Automatically) ##
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)); 

require(ggplot2)
require(cowplot)

survey22 <- readRDS("data_survey22_v7.rds")
icpp09 <- readRDS("data_icpp09_v7.rds")
icpp13 <- readRDS("data_icpp13_v7.rds")
wvs7 <- readRDS("wvs7_japan_v1.rds")

##########################
## Main DVs (Figure A8) ##
##########################

p09_1 <- ggplot(icpp09, aes(x=immigincrease)) + 
  geom_bar() + 
  labs(x=paste0("Admit More Foreigners", 
                "\n(Mean = ", round(mean(icpp09$immigincrease, na.rm=TRUE),3), ", ",
                "SD = ", round(sd(icpp09$immigincrease, na.rm=TRUE),3), ")"),  
       y="Frequency", title = "ICPP 2009") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5))
p09_1

p09_2 <- ggplot(icpp09, aes(x=foreignsuff)) + 
  geom_bar() + 
  labs(x=paste0("Integrate Foreign Residents", 
                "\n(Mean = ",round(mean(icpp09$foreignsuff, na.rm=TRUE),3), ", ",
                "SD = ", round(sd(icpp09$foreignsuff, na.rm=TRUE),3), ")"),  
       y="Frequency", title = "ICPP 2009") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5))
p09_2

p22_1 <- ggplot(survey22, aes(x=immigincrease)) + 
  geom_bar() + 
  labs(x=paste0("Admit More Foreigners", 
                "\n(Mean = ", round(mean(survey22$immigincrease, na.rm=TRUE),3), ", ",
                "SD = ", round(sd(survey22$immigincrease, na.rm=TRUE),3), ")"),  
       y="Frequency", title = "2022 Survey") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5))
p22_1

p22_2 <- ggplot(survey22, aes(x=foreignsuff)) + 
  geom_bar() + 
  labs(x=paste0("Integrate Foreign Residents", 
                "\n(Mean = ",round(mean(survey22$foreignsuff, na.rm=TRUE),3), ", ",
                "SD = ", round(sd(survey22$foreignsuff, na.rm=TRUE),3), ")"),  
       y="Frequency", title = "2022 Survey") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5))
p22_2

plot_grid(p09_1,p09_2,p22_1,p22_2, nrow=2)

ggsave("descr_dv.pdf", width=8, height=6)

#####################################################
## Main DVs, Mean and SD by age cohort (FIgure A9) ##
#####################################################

dvbyage <- rbind(do.call("rbind",tapply(icpp09$immigincrease, icpp09$cohort, 
                             function(k) c(mean(k,na.rm=T),sd(k,na.rm=TRUE)))),
      do.call("rbind",tapply(icpp09$foreignsuff, icpp09$cohort, 
                             function(k) c(mean(k,na.rm=T),sd(k,na.rm=TRUE)))),
      do.call("rbind",tapply(survey22$immigincrease, survey22$cohort, 
                             function(k) c(mean(k,na.rm=T),sd(k,na.rm=TRUE)))),
      do.call("rbind",tapply(survey22$foreignsuff, survey22$cohort, 
                             function(k) c(mean(k,na.rm=T),sd(k,na.rm=TRUE)))))
dvbyage
dvbyage <- as.data.frame(dvbyage)
colnames(dvbyage) <- c("mean","sd")
dvbyage$cohort <- factor(c(rep(levels(icpp09$cohort),2),rep(levels(survey22$cohort),2)), 
                         levels=rev(c(levels(survey22$cohort)[1],
                                      levels(icpp09$cohort),levels(survey22$cohort)[4])))
dvbyage$data <- factor(rep(c("ICPP 2009", "2022 Web Survey"), each=8),
                       levels=c("ICPP 2009", "2022 Web Survey"))
dvbyage$dv <- factor(rep(c("Admit More Foreigners", "Integrate Foreign Residents"), each=4),
                     levels=c("Admit More Foreigners", "Integrate Foreign Residents"))

ggplot(dvbyage, aes(x=cohort, y=mean, ymin=mean-sd, ymax=mean+sd)) + 
  geom_point() + 
  geom_errorbar(width=0.2) + 
  facet_grid(data~dv, scales = "free_y") + 
  coord_flip(ylim=c(0,1)) + 
  labs(y= "Mean ± One Standard Deviation", x=NULL) + 
  theme_bw()

ggsave("descr_dv_bycohort.pdf", width=7, height=5)

##################################
## Alternative DVs (Figure A10) ##
##################################

p09_1x <- ggplot(icpp09, aes(x=immigincrease_alt)) + 
  geom_histogram(bins=8, color="white") + 
  labs(x=paste0("Admit More Foreigners (Factor Score)", 
                "\n(Mean = ", round(mean(icpp09$immigincrease_alt, na.rm=TRUE),3), ", ",
                "SD = ", round(sd(icpp09$immigincrease_alt, na.rm=TRUE),3), ")"),  
       y="Frequency", title = "ICPP 2009") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5))
p09_1x

p09_2x <- ggplot(icpp09, aes(x=foreignrights)) + 
  geom_histogram(bins=8, color="white") + 
  labs(x=paste0("Integrate Foreign Residents (Factor Score)", 
                "\n(Mean = ",round(mean(icpp09$foreignrights, na.rm=TRUE),3), ", ",
                "SD = ", round(sd(icpp09$foreignrights, na.rm=TRUE),3), ")"),  
       y="Frequency", title = "ICPP 2009") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5))
p09_2x

p22_1x <- ggplot(survey22, aes(x=immigincrease_alt)) + 
  geom_bar() + 
  labs(x=paste0("Admit More Foreigners (Alternative)", 
                "\n(Mean = ", round(mean(survey22$immigincrease_alt, na.rm=TRUE),3), ", ",
                "SD = ", round(sd(survey22$immigincrease_alt, na.rm=TRUE),3), ")"),  
       y="Frequency", title = "2022 Survey") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5))
p22_1x

p22_2x <- ggplot(survey22, aes(x=foreignrights)) + 
  geom_histogram(bins=8, color="white") + 
  labs(x=paste0("Integrate Foreign Residents (Factor Score)", 
                "\n(Mean = ",round(mean(survey22$foreignrights, na.rm=TRUE),3), ", ",
                "SD = ", round(sd(survey22$foreignrights, na.rm=TRUE),3), ")"),  
       y="Frequency", title = "2022 Survey") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5))
p22_2x

p13_1x <- ggplot(icpp13, aes(x=immigincrease_alt)) + 
  geom_histogram(bins = 8, color="white") + 
  labs(x=paste0("Admit More Foreigners (Factor Score)", 
                "\n(Mean = ", round(mean(icpp13$immigincrease_alt, na.rm=TRUE),3), ", ",
                "SD = ", round(sd(icpp13$immigincrease_alt, na.rm=TRUE),3), ")"),  
       y="Frequency", title = "ICPP 2013") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5))
p13_1x

p13_2x <- ggplot(icpp13, aes(x=foreignrights)) + 
  geom_histogram(bins = 8, color="white") + 
  labs(x=paste0("Integrate Foreign Residents (Factor Score)", 
                "\n(Mean = ",round(mean(icpp13$foreignrights, na.rm=TRUE),3), ", ",
                "SD = ", round(sd(icpp13$foreignrights, na.rm=TRUE),3), ")"),  
       y="Frequency", title = "ICPP 2013") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5))
p13_2x

plot_grid(p09_1x,p09_2x,p22_1x,p22_2x,p13_1x,p13_2x, nrow=3)

ggsave("descr_dv_alt.pdf", width=8, height=9)


#############################
## WVS Wave 7 (Figure A11) ##
#############################

pwvs7_1 <- ggplot(wvs7, aes(x=immigincrease)) + 
  geom_bar() + 
  labs(x=paste0("Admit Foreign Workers", 
                "\n(Mean = ", round(mean(wvs7$immigincrease, na.rm=TRUE),3), ", ",
                "SD = ", round(sd(wvs7$immigincrease, na.rm=TRUE),3), ")"),  
       y="Frequency", title = "WVS Wave 7 (2019)") + 
  theme_bw() + 
  theme(plot.title = element_text(hjust=0.5))
pwvs7_1

ggsave("descr_dv_wvs7.pdf", width=5, height=4)
